#include <stdio.h>
#include <stdlib.h>
#include <time.h>

/**
 * 求解超定方程组组：
 *  2x + 3y = 5  （1） 
 *  x - 2y = 3    (2)  
 *  3x + y = 4    (3)
 *  
 *  转化为求使得误差:
 *       L(x, y) = 1/(2 * 3) ((2x + 3y - 5)^2 + (x - 2y -3)^2 + (3x + y -4)^2)
 *  最小化的优化问题。
 *
 *  梯度：▽L_(x, y) = 1/3 *(14x +7y - 25, 7x + 14y - 13)
 */

typedef double(*grad_fun)(double, double);
typedef double(*loss_fun)(double, double);

 double Loss(double x, double y)
 {
    return 1/6.0* (
        (2* x + 3* y - 5)* (2* x + 3* y - 5)
        + (x - 2 * y -3) * (x - 2 * y -3) 
        + (3* x + y -4) * (3* x + y -4));
 }

 double Lx(double x, double y)
 {
    return 1/3.0 * (14 * x + 7 * y - 25);
 }

 double Ly(double x, double y) 
 {
    return 1 / 3.0 * ( 7 * x + 14 * y - 13);
 }

 /**
  * @param x: x的初始值
  * @param y: y的初始值
  * @param partial_x: 计算误差关于x的偏导数的函数
  * @param partial_y: 计算误差关于y的偏导数的函数
  * @param loss_fn: 计算误差的函数
  * @param ox: 最终x的值
  * @param oy：最终y的值
  * @param lr: 学习率，控制沿梯度负方向前进的步幅
  * @param epochs: 迭代次数
  * 
  * @return 优化结束时的误差
  */
 double gradient_descent(
    double x, double y, 
    grad_fun partial_x, grad_fun partial_y, loss_fun loss_fn,
    double *ox, double *oy, 
    double lr, int epochs)
 {
    for (int i = 0; i < epochs; i++) {
        x += - lr * partial_x(x, y);
        y += - lr * partial_y(x, y);
        double loss = loss_fn(x, y);
        printf("epoch %d, loss = %.4f\n", i + 1, loss);
    }
    *ox = x;
    *oy = y;
    return loss_fn(x, y);
 }


 int main()
 {
    // 随机初始化x, y
	srand(time(NULL));
    double x = rand() % 100 * 1e-5;
    double y = rand() % 100 * 1e-5;

    printf("初始值 x = %.4f, y = %.4f, MSE = %.7f\n", x, y, Loss(x, y));
    double ret = gradient_descent(x, y, Lx, Ly, Loss, &x, &y, 1e-2, 500);
    printf("优化结果： x = %.4f, y = %.4f, MSE = %.7f\n", x, y, Loss(x, y));

    return 0;
 }
 